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Abstract In this paper, we present a numerical technique for performing Lie 
advection of arbitrary differential forms. Leveraging advances in high-resolution 
finite volume methods for scalar hyperbolic conservation laws, we first discretize 
the interior product (also called contraction) through integrals over Eulerian ap- 
proximations of extrusions. This, along with Cartan's homotopy formula and a 
discrete exterior derivative, can then be used to derive a discrete Lie derivative. 
The usefulness of this operator is demonstrated through the numerical advection of 
scalar fields and 1 -forms on regular grids. 
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1 Introduction 

Deeply-rooted assumptions about smoothness and differentiability of most continu- 
ous laws of mechanics often clash with the inherently discrete nature of computing 
on modem architectures. To overcome this difficulty, a vast number of compu- 
tational techniques have been proposed to discretize differential equations, and 
numerical analysis is used to prove properties such as stability, accuracy, and con- 
vergence. However, many key properties of a mechanical system are characterized 
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by its symmetries and invariants (e.g., momenta), and preserving these features 
in the computational realm can be of paramount importance 1 19 1, independent of 
the order of accuracy used in the computations. To this end, geometrically-derived 
techniques have recently emerged as valuable alternatives to traditional, purely 
numerical-analytic approaches. In particular, the use of differential forms and their 
discretization as cochains has been advocated in a number of applications such as 
electromagnetism 16114111231 . discrete mechanics (STJ, and even fluids |[T4l . 

In this paper we introduce a finite volume based technique for solving the discrete 
Lie advection equation, ubiquitous in most advection phenomena: 

--^Cxu; = 0, (1) 
at 

where is an arbitrary discrete differential /c-form ISllSlfTTl defined on a discrete 
manifold, and X is a discrete vector field living on this manifold. Our numerical 
approach stems from the observation, developed in this paper, that the computa- 
tional treatment of discrete differential forms share striking similarities with finite 
volume techniques |27 | and scalar advection techniques used in level sets I37II36I . 
Consequently, we present a discrete interior product (or contraction) computed us- 
ing any of the /c-dimensional finite volume methods readily available, from which 
we derive a numerical approximation of the spatial Lie derivative Cx using a 
combinatorial exterior derivative. 

1.1 Background on the Lie Derivative 

The notion of Lie derivative Cx in Elie Cartan's Exterior Calculus |10 | extends 
the usual concept of the derivative of a function along a vector field X. Although a 
formal definition of this operator can be made purely algebraically (see 1 1 1, §5.3), 
its nature is better elucidated from a dynamical perspective [1] (§5.4). Consequently, 
the spatial Lie derivative (along with its closely related time-dependent version) 
is an underlying element in all areas of mechanics: for example, the rate of strain 
tensor in elasticity and the vorticity advection equation in fluid dynamics are both 
nicely described using Lie derivatives. 

A common context where a Lie derivative is used to describe a physical evolution 
is in the advection of scalar fields: a scalar field p being advected in a vector field 
V can be written as: dp/dt -\- C^p = 0. The case of divergence-free vector 
fields (i.e., V • V = 0) has been the subject of extensive investigation over the past 
several decades leading to several numerical schemes for solving these types of 
hyperbolic conservation laws in various applications (see, e.g., II38[I12[|42II25[I22[ 
[T4l[T3]| ). Chief among them are the so-called finite volume methods ||27]| . including 
upwind, ENO, WENO, and high-resolution techniques. Unlike finite difference 
techniques based on point values (e.g., I 13 Bq1|291| ), such methods often resort to 
the conservative form of the advection equation and rely on cell averages and the 
integrated fluxes in between. The integral nature of these finite volume techniques 
will be particularly suitable in our context, as it matches the foundations behind 
discrete versions of exterior calculus EO. 
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While finite volume schemes have been successfully used for over a decade, they 
have been used almost solely to advect scalar fields, be they functions or densities, 
or systems thereof (e.g., components of tensor fields). To the authors' knowledge. 
Lie advection of inherently non-scalar entities such as differential forms has yet to 
benefit from these advances, as differential forms are not Lie advected in the same 
manner as scalar fields. 

1.2 Emergence of Structure-Preserving Computations 

Concurrent to the development of high-resolution methods for scalar advection, 
structure-preserving geometric computational methods have emerged, gaining ac- 
ceptance among engineers as well as mathematicians |2|. Computational electro- 
magnetism 1 6, 41 1, mimetic (or natural) discretizations |35,5|, and more recently 
Discrete Exterior Calculus (DEC, I24II11I ) and Finite Element Exterior Calculus 
(FEEC, 1 3 1) have all proposed similar discrete structures that discretely preserve 
vector calculus identities to obtain improved numerics. In particular, the relevance 
of exterior calculus (Cartan's calculus of differential forms ifTOl ) and algebraic 
topology (see, for instance, 1 33 1) to computations came to light. 

Exterior calculus is a concise formalism to express differential and integral equa- 
tions on smooth and curved spaces in a consistent manner, while revealing the 
geometrical invariants at play. At its root is the notion of differential forms, denot- 
ing antisymmetric tensors of arbitrary order. As integration of differential forms 
is an abstraction of the measurement process, this calculus of forms provides an 
intrinsic, coordinate-free approach particularly relevant to concisely describe a 
multitude of physical models that make heavy use of line, surface and volume 
integrals ISlfTl BOII 1 6l[32l l9l[T7l . Similarly, many physical measurements, such as 
fluxes, are performed as specific local integrations over a small surface of the mea- 
suring instrument. Pointwise evaluation of such quantities does not have physical 
meaning; instead, one should manipulate those quantities only as geometrically- 
meaningful entities integrated over appropriate submanifolds — these entities and 
their geometric properties are embodied in discrete differential forms. 

Algebraic topology, specifically the notion of chains and cochains (see, e.g., 1441 
[33 1, has been used to provide a natural discretization of these differential forms 
and to emulate exterior calculus on finite grids: a set of values on vertices, edges, 
faces, and cells are proper discrete versions of respectively pointwise functions, 
line integrals, surface integrals, and volume integrals. This point of view is entirely 
compatible with the treatment of volume integrals in finite volume methods, or 
scalar functions in finite element methods |5|; but it also involves the edge ele- 
ments and facet elements as introduced in E&M as special i^div and i^curi basis 
elements |34|. Equipped with such discrete forms of arbitrary degree. Stokes' the- 
orem connecting differentiation and integration is automatically enforced if one 
thinks of differentiation as the dual of the boundary operator — a particularly simple 
operator on meshes. With these basic building blocks, important structures and 
invariants of the continuous setting directly carry over to the discrete world, culmi- 
nating in a discrete Hodge theory (see recent progress in |4 1). As a consequence, 
such a discrete exterior calculus has, as we have mentioned, already proven use- 
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ful in many areas such as electromagnetism |[6ll4T1l . fluid simulation lfT4ll . surface 
parameterization Iil8i , and remeshing of surfaces B3ll to mention a few. 

Despite this previous work, the contraction and Lie derivative of arbitrary discrete 
forms — two important operators in exterior calculus — have received very little 
attention, with a few exceptions. The approach in |7 1 (that we will review in §3.1[ ) is 
to exploit the duality between the extrusion and contraction operators, resulting in 
an integral definition of the interior product that fits the existing foundations. While 
a discrete contraction was derived using linear "Whitney" elements, no method to 
achieve low numerical diffusion and/or high resolution was proposed. Furthermore, 
the Lie derivative was not discussed. More recently Heumann and Hiptmair 12111 
leveraged this work to suggest an approach similar to ours in a finite element 
framework for Lie advection of forms of arbitrary degree, however only 0-forms 
were analyzed. 

1.3 Contributions 

In this paper we extend the discrete exterior calculus machinery by introducing 
discretizations of contraction and Lie advection with low numerical diffusion. 
Our work can also be seen as an extension of classical numerical techniques for 
hyperbolic conservation laws to handle advection of arbitrary discrete differential 
forms. In particular, we will show that our scheme in 3D is a generalization of 
finite volume techniques where not only cell-averages are used, but Silso face- and 
edge-averages, as well as vertex values. 

2 Mathematical Tools 

Before introducing our contribution, we briefly review the existing mathematical 
tools we will need in order to derive a discrete Lie advection: after discussing our 
setup, we describe the necessary operators of Discrete Exterior Calculus, before 
briefly reviewing the foundations of finite volume methods for advection. In this 
paper continuous quantities and operators are distinguished from their discrete 
counterparts through a bold typeface. 

2. 1 Discrete Setup 

Space Discretization. Throughout the exposition of our approach, we assume a 
regular Cartesian grid discretization of space. This grid forms an orientable 3- 
manifold cell complex K = (F, F, C) with vertex set F = {vi}, edge set 
E = {cij}, as well as face set F and cell set C. Each cell, face and edge is assigned 
an arbitrary yet fixed intrinsic orientation, while vertices and cells always have a 
positive orientation. By convention, if a particular edge Cij is positively oriented 
then Cji refers to the same edge with negative orientation, and similar rules apply 
for higher dimensional mesh elements given even vs. odd permutations of their 
vertex indexing. 

Boundary Operators. Assuming that mesh elements in K are enumerated with an 
arbitrary (but fixed) indexing, the incidence matrices of K then define the boundary 
operators. For example, we let denote the \V\ x \E\ matrix with {d^)ve = 1 
(resp., —1) if vertex v is incident to edge e and the edge orientation points towards 
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(resp., away from) v, and zero otherwise. Similarly, denotes the |£^| x |F| 
incidence matrix of edges to faces with {d^)ef = 1 (resp., —1) if edge e is incident 
to face / and their orientations agree (resp., disagree), and zero otherwise. The 
incidence matrix of faces to cells is defined in a similar way. See 1331 for details. 

2.2 Calculus of Discrete Forms 

Guided by Cartan's exterior calculus of differential forms on smooth manifolds, 
DEC offers a calculus on discrete manifolds that maintains the covariant nature of 
the quantities involved. 

Chains and Cochains. At the core of this computational tool is the notion of chains, 
defined as a linear combination of mesh elements; a 0-chain is a weighted sum 
of vertices, a 1 -chain is a weighted sum of edges, etc. Since each /c-dimensional 
cell has a well-defined notion of boundary (in fact its boundary is a chain itself; 
the boundary of a face, for example, is the signed sum of its edges), the boundary 
operator naturally extends to chains by linearity. A discrete form is simply defined 
as the dual of a chain, or cochain, a linear mapping that assigns each chain a 
real number. Thus, a 0-cochain (that we will abusively call a 0-form sometimes) 
amounts to one value per 0-dimensional cell, such that any 0-chain can naturally 
pair with this cochain. More generally, /c-cochains are defined by one value per 
/c-cell, and they naturally pair with /c-chains. The resulting pairing of a /c-cochain 
and a /c-chain is the discrete equivalent of the integration of a continuous 
/c-form over a /c-dimensional submanifold tr/ei 



While attractive from a computational perspective due to their conceptual simplicity 
and elegance, the chain and cochain representations are also deeply rooted in a 
theoretical framework defined by H. Whitney |44|, who introduced the Whitney 
and deRham maps that establish an isomorphism between the cohomology of 
simplicial cochains and the cohomology of Lipschitz differential forms. With these 
theoretical foundations, chains and cochains are used as basic building blocks 
for direct discretizations of important geometric structures such as the deRham 
complex through the introduction of two simple operators. 

Discrete Exterior Derivative. The differential d (called exterior derivative) is an 
existing exterior calculus operator that we will need in our construction of a Lie 
derivative. The discrete derivative d is constructed to satisfy Stokes' theorem, which 
elucidates the duality between the exterior derivative and the boundary operator. In 
the continuous sense, it is written 



Consequently, if a is a discrete differential /c-form, then the (/c+l)-form da is 
defined on any (/c+l)-chain a by 





(2) 



(da, a) = (a, da) , 



(3) 
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where da is the (/c-chain) boundary of a, as defined in ^2.1 Thus the discrete dif- 
ferential d, mapping /c-forms to (/c+l)-forms, is given by the co-boundary operator, 
the transpose of the signed incidence matrices of the complex K\ do = {d^)^ 
maps 0-forms to 1-forms, di = (9^)^ maps 1-forms to 2-forms, and more gen- 
erally in nD, dk = (9^+^)^. In relation to standard 3D vector calculus, this can 
be seen as do = V,di = Vx, and d2 = The fact that the boundary of 
a boundary is empty results in dd = 0, which in turn corresponds to the vector 
calculus facts that VxV = V- Vx = 0. Notice that this operator is defined 
purely combinatorially, and thus does not need a high-order definition, unlike the 
operators we will introduce later. 

2.3 Principles of Finite Volumes 

Given the integral representation of discrete forms used in the previous section, a 
last numerical tool we will need is a method for computing solutions to advection 
problems in integral form. Finite volume methods were developed for exactly this 
purpose, and while we now provide a brief overview of this general procedure for 
completeness, we refer the reader to 1 27 1 and references therein for further details 
and applications. One approach of finite volume schemes is to advect a function 
u{x) by a velocity field v{x) using a Reconstruct-Evolve- Average (REA) approach. 
In one dimension, we can define the cell average of a function u{x) over cell Ci 
with width Z\x as 



u{x) dx i = 1, 2, . . . , A^. 



Given k adjacent cell averages, the method will reconstruct a function such that the 
average ofp{x) in each of the k cells is equal to the average ofu{x) in those cells. 
High-resolution methods attempt to build a reconstruction such that it has only high- 
order error terms in smooth regions, while lowering the order of the reconstruction 
in favor of avoiding oscillations in regions with discontinuities like shocks. Such 
adaptation can be done through the use of slope limiters or by changing stencil sizes 
using essentially non-oscillatory (ENO) and related methods. This reconstruction 
can then be evolved by the velocity field and averaged back onto the Eulerian grid. 

Another variant of finite volume methods is one that computes fluxes through cell 
boundaries. Employing Stokes' theorem, the REA approach can be implemented 
by computing only the integral of the reconstruction which is evolved through each 
face, and then differencing the incoming and outgoing integrated fluxes of each 
cell to determine its net change in density. It is this flux differencing approach 
that will be most convenient for deriving our discrete contraction operator, due 
to the observation that the net flux through a face induced by evolving a function 
forward in a velocity field is equal to the flux through the face induced by evolving 
the face backwards through the same velocity field. This second interpretation 
of the integrated flux is the same as computing the integral of the function over 
an extrusion of the face in the velocity field, as will be seen in the next section, 
and therefore we may use any of the wide range of finite volume methods to 
approximate integrals over extruded faces. 
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3 Discrete Interior Product and Discrete Lie Derivative 

In keeping with the foundations of Discrete Exterior Calculus, we present the 
continuous interior product and Lie derivative operators in their "integral" form, i.e., 
we present continuous definitions of ix^ and Cx(-^ integrated over infinitesimal 
submanifolds: these integral forms will be particularly amenable to discretization 
via finite volume methods and DEC as we discussed earlier. 

3.1 Towards a Dynamic Definition of Lie Derivative 

Interior Product through Extrusion. As pointed out by |7 1, the extrusion of objects 
under the flow of a vector field can be used to give an intuitive dynamic definition of 
the interior product. If Al is an n-dimensional smooth manifold and X G X( Al) a 
smooth (tangent) vector field on the manifold, let <S be a /c-dimensional submanifold 
on M with k < n. The flow if of the vector field X is simply a function (p : 
M X R ^ M consistent with the one-parameter (time) group structure, that is, 
such that (p{(p{S, t),s) = (p{S, s ^ t) with (p{S, 0) = S for all s,t eR. Now 
imagine that S is carried by this flow of X for a time t; we denote the resultant 
"flowed-out" submanifold Sx {t), which is equivalent to the image of S under the 
mapping (p, i.e., Sx{t) = cp{S^t). The extrusion Ex{S^t) is then the (k-\-l)- 
dimensional submanifold formed by the advection of S over the time t to its final 
position Sx {t) : it is the "extruded" (or "swept out") submanifold. This can be 
expressed formally as a union of flowed-out manifolds, 

Ex{S,t)= y Sx{r) 

re[o,t] 

where the orientation of Ex (<5, t) is defined such that 

dEx=Sx{t)-S-Ex{dS,t). 

These geometric notions are visualized in Figure [T] where the submanifold S is 
presented as a 1-dimensional curve, flowed out to form Sx{t), or alternatively, 
extruded to form Ex (<5, 

Using this setup, the interior product ix of a time-independent form u; evaluated 
on S can now be defined through one of its most crucial properties, i.e., as the 
instantaneous change of evaluated on Ex {^jt), or more formally. 



uj. (4) 

While this equation is coherent with the discrete spatial picture, for the discrete Lie 
advection we will also wish to integrate ix^ over a small time interval. Hence, 
by taking the integral of both sides of Eq. ^ over the interval [0, At], the first 
fundamental theorem of calculus gives us 



I I ixuj dt= [ uj, (5) 

JO JSxit) J Ex (S, At) 



which will be used later on for the discretization of the time-integrated interior 
product. 
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X 

Ex{dS,±t) 
dS 



Fig. 1 Geometric interpretation of the Lie derivative Cx(^ of a differential form in 
the direction of vector field X: (a) for a backwards advection in time of an edge S (referred 
to as upwind extrusion), and (b) for a forward advection of S. Notice the orientation of the 
two extrusions are opposite, and depend on the direction of the velocity field. 

Algebraic and Flowed-out Lie Derivative. Using a similar setup, we can formulate 
a definition of Lie derivative based on the flowed-out submanifold Sx{t). Re- 
member that the Lie derivative is a generalization of the directional derivative to 
tensors, intuitively describing the change of u: in the direction of X. In fact, the 
Lie derivative Cx(-^ evaluated on S is equivalent to the instantaneous change of u; 
evaluated on Sx {t), formally expressed by 



Cx<^ 



d 
dt 



(6) 



t=0 -^Sxit 



as a direct consequence of the Lie derivative theorem |[Q(Theorem 6.4.1). As 
before, we can integrate Eq. ([6]) over a small time interval [0, At], applying the 
Newton-Leibnitz formula to find 



r I 

Jo Js 



Cxu^ 



dt= uj- i 



(7) 



Note that the formulation above, discretized using a semi-Lagrangian method, has 
been used, e.g., by 1 14| to advect fluid vorticity; in that case the right hand side 
of Eq. ^ was evaluated by looking at the circulation through the boundary of 
the "backtracked" manifold. Rather than following their approach, we revert to 
discretizing the dynamic definition of the interior product in Eq. ^ instead, and 
later constructing the Lie derivative algebraically. The primary motivation behind 
this modification is one of effective numerical implementation: we can apply a 
dimension-by-dimension finite volume scheme to obtain an approximation of the 
interior product, while the alternative — computing integrals of approximated uj over 
Sx (t) as required by a discrete version of Eq. ^ — is comparatively cumbersome. 
Also, by building on top of standard finite volume schemes the solvers can leverage 
pre-existing code, such as CLAWR\CK (261, without requiring modification. 

We now show how the Lie derivative and the interior product are linked through a 
simple algebraic relation known as Cartan's homotopy formula. In particular, this 
derivation (using Figure [T] as a reference) requires repeated application of Stokes' 
theorem from Eq. ([2]). 



lim ^ 

At^O At 



1 r^^ r 

^tJo Js 



Cx(^ 

ISx{t) 



dt = lim — - 

At^O At 



JSx{At) Js 



(8) 
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lim —— 

At^O At 



J Ex (S, At) Je 



UJ 

Ex {OS, At) 



= / ixdu?+ / ] 
Js JdS 



Js Js 



(9) 
(10) 
(11) 



The submanifolds S and Sx {^i) form a portion of the boundary of Ex (<5, ^t)- 
Therefore, by Stokes', we can evaluate duj on the extrusion and subtract off the 
other portions of dEx (<5, ^t) to obtain the desired quantity. This is how we pro- 
ceed from Eq. ^ to Eq. ^ of the proof. The following line, Eq. ([TO]), is obtained 
by applying the dynamic definition of the interior product given in Eq. ([5]) to each 
of two terms, leading us to our final result in Eq. ( pTj ) through one final applica- 
tion of Stokes' theorem. What we have obtained is the Lie derivative expressed 
algebraically in terms of the exterior derivative and interior product. Notice that 
Eq. ( pT] ) is the integral form of the celebrated identity called Cartan's homotopy (or 
magic) formula, most frequently written as 



Cx^ = ixdu? + dix<^. 



(12) 



By defining our discrete Lie derivative through this relation, we ensure the algebraic 
definition holds true in the discrete sense by construction. It also implies that the Lie 
derivative can be directly defined through interior product and exterior derivative, 
without the need for its own discrete definition. 

Upwinding the Extrusion We may rewrite the above notions using an "upwinded" 
extrusion {i.e., a cell extruded backwards in time) as well (see Fig.[T|i). For example, 
Eq. ^ can be rewritten as 



s 



(13) 



t=0 JEx{S-t) 



While this does not change the instantaneous value of the contraction, integrating 
Eq. ([TSj over the time interval [0, At] now gives us 



/ / ixuj dt = - [ 

Jo Jsxit) Je 



UJ. 



(14) 



lEx{S-At) 



Similar treatment for the remainder of the above can be done and Cartan's formula 
can be derived the same way, however by using these definitions in our following 
discretization we will obtain computations over upwinded regions equivalent to 
those computed by finite volume methods. 

3.2 Discrete Interior Product 

A discrete interior product is computed by exploiting the principles of Eq. ([5j and 
applying the finite volume machinery. Given a discrete k-ioxm a and a discrete 
vector field X, the interior product is approximated by extruding backwards in 



10 



R Mullen et al. 



time every (A:-l)-dimensional cell a of the computational domain to form a new k- 
dimensional cell Ex{cr-, — At). Evaluating the integral of a over the extrusion and 
assigning the resulting value to the original cell cr yields the mapping (ix<^, cr) in- 
tegrated over a time step At. This procedure, once applied to all (/c-l)-dimensional 
cells, gives the desired discrete (/c-l)-form \xOi- 



Fig. 2 Approximating Extrusions: In the discrete setting, the extrusion of a (/c-l)- 
dimensional manifold (A:=l on left, 2 on right) is approximated by projecting the Lagrangian 



K-dimensional Splitting. One option for computing this integral would be to do 
an n-dimensional reconstruction of a, perform a Lagrangian advection of the cell 
cr to determine Ex{cf-, — At), and then algebraically or numerically computing 
the integral of the reconstructed a over this extrusion. In fact, this is the idea 
behind the approach suggested in ||2T1 . However, with the exception of when /c = n, 
such an approach does not allow us to directly leverage finite volume methods, 
as performing an n-dimensional reconstruction of a form given only integrals 
over /c-dimensional submanifolds would require a more general finite element 
framework. For simplicity and ease of implementation we avoid this generalization 
and instead resort to projecting the extrusion onto the grid-aligned /c-dimensional 
subspaces and then applying a /c-dimensional finite volume method to each of the 
(^) projections. The integrals over the extrusion of cr from each dimension are then 
summed. Again, note that in the special case of /c = n no projection is required and 
we are left exactly with an n-dimensional finite volume scheme. We have found 
that this splitting combined with a high-resolution finite volume method, despite 
imposing at most first order accuracy, can still give high quality results with low 
numerical diffusion, while being able to leverage existing finite volume solvers 
without modification. However, if truly higher order is required then a full-blown 
finite element method would most likely be required | 21 1. 

Finite Volume Evaluation. As hinted at in Section |2.3| we notice that the time 
integral of the flux of a density field being advected through a submanifold cr is 
equivalent to the integral of the density field over the backwards extrusion of cr 
over the same amount of time. In fact, some finite volume methods are derived 
using this interpretation, doing a reconstruction of the density field, approximating 
the extrusion, and integrating the reconstruction over this. However, many others 
are explained by computing a numerical flux per face, and then multiplying by 
the time step At: this is still an approximation of the integral over the extrusion, 
taking the reconstruction to be a constant (the numerical flux divided by v) and 
the extrusion having length vAt. Indeed the right hand side of Eq. ^ can be seen 
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as analogous to the numerical flux, after which Eq. ^ becomes the relationship 
between integrating the flux over time and the form over the extrusion. Hence we 
may use any of the finite volume methods for /c-dimensional density advection 
problems when computing the contraction of a /c-form. The only difference here 
is that rather than applying Stokes theorem and summing the contributions back 
to the original /c-cell (which will be done by the discrete exterior derivative in the 
dix^ term of the Lie derivative), the contraction simply stores the values on the 
(/c-l)-cells, without the final sum. 

3.3 Discrete Lie Advection 

We now have all the ingredients to introduce a discrete Lie advection. Given a k- 
form a, we compute the (/c+l)-fornida by applying the transpose of the incidence 



matrix 9^+^ to a as detailed in {22 We then compute the /c-form ix(da), and 
the (/c-l)-form By applying d to the latter form and summing the resulting 
/c-form with the other interior product, we finally get an approximation of Cartan's 
homotopy formula of the Lie derivative. An explicit example of this will be given 
in the next section to better illustrate the process and details. 

4 Applications and Results 

We now present a few direct applications of this discrete Lie advection scheme. In 
our tests we used upwinding one-dimensional WENO schemes for our contraction 
operator, splitting even the /c-dimensional problems into multiple one-dimensional 
ones. We found that when using high-resolution WENO schemes we could obtain 
quality results with little numerical smearing despite this dimensional splitting. 

A Note on Vector Fields. In this section we assume that vector fields are discretized 
by storing their flux {i.e., contraction with the volume form) on all the (n — 1)- 
dimensional cells of a nD regular grid, much like the Marker- And-Cell "staggered" 
grid setup 120]. Evaluation of the vector fields at lower dimensional cells is done 
through simple averaging of adjacent discrete fluxes. We pick this setup as it is one 
of the most commonly-used representations, but the vector fields can be given in 
arbitrary form with only minor implementation changes. 

4.1 Volume Forms and 0-Forms 

Applying our approach to volume forms (n-forms in n dimensions) we have 

Note that ix^ is the numerical flux computed by the chosen n-dimensional finite 
volume scheme while d will then just assign the appropriate sign of this flux to each 
cell's update, and hence we trivially arrive at the chosen finite volume scheme with 
no modification. Similarly, applying this approach to 0-forms results in well-known 
finite difference advection schemes of scalar fields. Indeed, we have in this case 

Cx^ = ixdct? + dix^u? = ixdcj 



as the contraction of a 0-form vanishes. We are thus left with du; computing 
standard finite differences of a node-based scalar field on edges, and ix then doing 
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Fig. 3 Grid setup: Indexing and location of the various quantities stored on different parts 
of the grid. Arrows indicate the orientation of the edges. 

componentwise upwind integration of reconstructions of these derivatives. Such 
techniques are common in scalar field advection, for example in the advection of 
level sets, and we refer the reader to L36..37J and references therein for examples. 

4.2 Advecting a 1-Form in 2D 

The novelty of this approach comes when applied to /c-forms in n dimensions with 
n > /c > 0. We first demonstrate the simplest such application of our method 
by advecting a 1-form by a static velocity field in 2D using the simple piecewise- 
constant up winding finite volume advection. To illustrate the general approach we 
will explicitly write out the algorithm for this case. We will assume the velocity 
X is everywhere positive in both x and y components to simplify the upwinding, 
and and X^ will be used to represent the integrated flux through vertical and 
horizontal edges respectively. 

Suppose we have a regular two-dimensional grid with square cells of size K^, and 
with each horizontal edge oriented in the positive x direction and each vertical edge 
oriented in the positive y direction and numbered according to Fig [3] A discrete 
1-form u) is represented by its integral along each edge. Due to the Cartesian 
nature of the grid, this implies that the dx component of the form will be stored on 
horizontal edges and the dy component will be stored on vertical edges, and we 
represent these scalars as ^ and u)\^- for the integrals along the (z, j) horizontal 
and vertical edge respectively. The discrete exterior derivative integrated over cell 
(z, j), {du))i^j, consists of the signed sum of uj over cell {i^jYs boundary edges, 
namely 



y 



Using piecewise-constant upwind advection, and remembering the assumption 
of positivity of the components of X, we may now compute ixdu? over a time 



Discrete Lie Advection of Differential Forms 



13 




(a) (b) (c) 

Fig. 4 1-Form advection: (a) A piecewise-constant form (dy within a rectangular shape, 
outside) is advected in a constant velocity field (X = (1,1), blue arrow) on a unit square 
periodic domain with a grid resolution of 48^ and a time step dt = 10~^. (b) Because the 
domain is periodic, the form should be advected back to its original position after Is (1000 
steps); however, our numerical method with a piecewise constant upwind finite volume 
scheme results in considerable smearing instead, (c) Using a high-resolution scheme (here, 
WENO-7) as the basic component of our form advection procedure significantly reduces 
smearing artifacts (same number of steps and step size). 

interval At for the horizontal and vertical edges (z, j) as 



(ixdu;)f,^. = -^X^(du;K,_i 



(15) 



Note the sign difference is due to the orientation of the extrusions, and would be 
different if the velocity field changed sign (see Fig.[T]). To compute the second half 
of Cartan's formula we must now compute ix^ at nodes, and then difference them 
along the edges. Using dimension splitting, as well averaging the velocity field 
from edges to get values at nodes, we get for node (z, j) 

(ixu;)i.i = ^ {{X!^j + + {XI, + . (16) 

We may now trivially compute dix^ for edges as 



(dix<^> 



(ix<^) 



{d\x^)l. = (ixcc?)i,i+i - {ix(^)i,j- 
Cartan's formula and the definition of Lie advection now lead us to obtain 

rAt 



-f 

Jo 



Cx^dt 



discretized as the update rule 



-{ixduj)lj - (dix^)f,^- 
-(}xduj)lj - {d\x^)l.. 
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Fig. 5 High-order Advection: in a vortical vector field (left) typically used for scalar 
advection, a piecewise-constant form is advected on a unit square periodic domain with a 
grid resolution of 48^ and a time step dt = 10"^ for 0, 200, and 400 steps (top), 600, 800, 
and 1000 steps (bottom). 

A First Example. An example of this low-order scheme can be seen in Figure Qa-b) 
where we advect a piecewise constant 1-form by a constant vector field X = (1,1) 
in a periodic domain. Advecting the form forward in this velocity field for a time of 
Is brings the form back to its original position in the continuous case; however, our 
numerical scheme proves very diffusive, as expected on discontinuous forms. We 
can however measure the error of our scheme by comparison with initial conditions 
as a function of the grid resolution with appropriately scaled time step sizes. To 
measure the error we recall the norm of a /c-form cv is defined over a smooth 
manifold Al as 

\^\p = 

and (•, -j^vH is the scalar product of /c-forms defined by the Riemannian metric, and 
d/ii is its associated volume form. We hence define the 1- and 2-norms of discrete 
1 -forms on a 2D regular grid with spacing h as 

\u;U = hY,{K,\ + K,\) 

for simplicity, but we found using more sophisticated discretizations of the norms 
all yielded similar results. Figure [6jc) shows the error plot in Li and L2 norms 
of this simple example under power-of-two refinement, confirming the first-order 
accuracy of our approach. 



/ 



where 



Discrete Lie Advection of Differential Forms 



15 



High-Resolution Methods Note that had we chosen to leverage more sophisticated 
finite volume solvers in the previous example, the only changes would occur in 
Equations ( p3] ) and ([16]) which would use the new numerical flux for computing the 
discrete contraction: any 2D method could be used for Equation ( p3] ), while a ID 
method is required for Equation ([16]). Due to the dimensional splitting obtaining 
higher order schemes is not easy, but for many application the order of accuracy is 
not always the most important thing. In particular, in the presence of discontinuous 
solutions high-resolution methods are often preferred for their ability to better 
preserve discontinuities and reduce diffusion. To test the utility and effectiveness 
of such schemes applied to forms, we compare the piecewise-constant up winding 
method from the previous section with a Finite Volume T^'^-order ID WENO upwind 
scheme (see an overview of FV-WENO methods in [39]). Figure |4]^c) shows the 
high-resolution finite volume scheme does a much better job at preserving the 
discontinuities, despite both methods being of the same order of accuracy for this 
discontinuous initial form (Figure |6];c)). 





Grid Resolution (log scale) 

(c) 



Fig. 6 Error Plots and Convergence 
Rates: We provide error plots in Li and 
L2 norms for power-of-two refinements 
for (a) advection of a smooth form in a 
constant vector field, (b) advection of a 
smooth form in the vortical vector field 
of Figure [5jleft), and (c) advection of the 
discontinuous form used in Figure [5] The 
black bold segments indicate a slope of 1. 



Accuracy To further demonstrate that properties from the underlying finite volume 
schemes chosen (including their accuracy) carry over to the advection of forms, we 
provide additional numerical tests. In Figure [5] we advected a simple discontinuous 
1-form in a vortical shearing vector field (Rudman vortex, left) on a 48x48 grid 
representing a periodic domain. As expected, the form is advected in a spiral-like 
fashion. By advecting the shape back in the the negated velocity field for the same 
amount of time, we can derive error plots to compare the Li and L2 norms for 
this example under refinement of the grid; see Figure |6]for other convergence tests 
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when a first-order upwind method or a WENO-7 method is used in our numerical 
technique. 

4.3 Properties 

It is easy to show that our discrete Lie derivative will commute with the discrete 
exterior derivative as in the continuous case, by using Cartan's formula and a 
discrete exterior derivative which satisfies dd = since we have 

dCx^ = d(ixdu; + dix^) 
= dixdcj 
= (dix + ixd)du; 
= Z^xdu?. 

This commutativity does not depend on any properties of the discrete contraction 
and therefore holds regardless of the underlying finite volume scheme chosen. A 
useful consequence of this fact is that the discrete Lie advection of closed forms 
will remain closed by construction; i.e., the advection of a gradient (resp., curl) 
field will remain a gradient (resp., curl) field. 

Unfortunately other properties of the Lie derivative do not carry over to the discrete 
picture as easily. The product rule for wedge products, for example, does not 
hold for the discrete wedge products defined in 1241 . although perhaps a different 
discretization of the wedge product may prove otherwise. However, the nonlinearity 
of the discrete contraction operator along with the upwinding potentially picking 
different directions on simplices and their subsimplices makes designing a discrete 
analog satisfying this continuous property challenging. 

5 Conclusions 

In this paper we have introduced an extension of classical finite volume techniques 
for hyperbolic conservation laws to handle arbitrary discrete forms. A class of first- 
order finite- volume-based discretizations of both contraction and Lie derivative 
of arbitrary forms was presented, extending Discrete Exterior Calculus to include 
approximations to these operators. Low numerical diffusion is attainable through 
the use of high-resolution finite volume methods. The advection of forms and vector 
fields are applicable in a multitude of problems, including conservative interface 
advection and conservative vorticity evolution. 

Although finite volume methods can offer high resolution at a relatively low com- 
putational cost, numerical diffusion is still present and can accumulate over time. 
In addition, the numerical scheme we presented is not variational in nature, i.e., 
it is not (a priori) derived from a variational principle. These limitations are good 
motivations for future work. 

While we have given numerical evidence demonstrating the manner in which 
resolution and accuracy is inherited from the underlying finite volume scheme, a 
formal analysis of stability and convergence remains to be performed. In particular, 
it is desirable to understand what the stability of the underlying one-dimensional 
scheme implies about the resulting stability of our method. Although norms for 
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discrete differential forms have been defined, such tools are not always suitable for 
the analysis of nonlinear methods in multiple dimensions, even for scalar advection. 
Such an investigation would thus be an important next step for the present work. 

In the future, we also expect that extensions can be made to make truly high-order 
and high-resolution discretizations of the contraction and Lie derivative through 
n-dimensional reconstructions of /c-forms and extrusions. In particular, this would 
greatly facilitate the extension to simplicial meshes. While there has been recent 
progress on high-order schemes for triangular meshes |45 , 28], these are not directly 
applicable for contractions of arbitrary forms. Hence, despite having a well-defined 
discrete exterior derivative for simplicial meshes, such an extension will require 
more than just 0- or n-form advection schemes, as our current dimension splitting 
approach to generalize such schemes for contraction of forms of other degrees does 
not immediately apply on non-rectangular meshes. 
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